Bagnold scaling, density plateau, and kinetic theory analysis of dense granular flow 
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We investigate the bulk rheology of dense granular flow down a rough slope, where the density 
profile has been found to show a plateau except for the boundary layers in simulations [Silbert et 
al, Phys. Rev. E 64, 05f302 (2001)]. It is demonstrated that both the Bagnold scaling and the 
framework of kinetic theory are applicable in the bulk, which allows us to extract the constitutive 
relations from simulation data. The detailed comparison of our data with the kinetic theory shows 
quantitative agreement for the normal and shear stresses, but there exists slight discrepancy in the 
energy dissipation, which causes rather large disagreement in the kinetic theory analysis of the flow. 
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Flowing granular material behaves like a fluid, but 
comprehensive understanding of its rheology is still far 
from complete. In the low-density regime with large 
shear rate, grains interact through instantaneous colli- 
sions and are described by a hydrodynamic model based 
on kinetic theory of inelastic hard spheres . As the sys- 
tem becomes denser, the independent collision assump- 
tion becomes questionable and the one particle distribu- 
tion of grain velocity may not be characterized by a small 
number of parameters or temperatures. When grains are 
nearly closed packed, they may experience enduring con- 
tacts, and the system behaves as in plastic deformation. 

One of a few established laws that hold for granular 
rheology up to the relatively dense regime is the Bagnold 
scaling , which states that the shear stress is propor- 
tional to the square of the strain rate. In fact, this is 
the only possible form for the stress in the flow of rigid 
grains characterized by the shear rate 7 and the packing 
fraction v, because 7 -1 is the only time scale. Simple 
dimensional consideration gives the Bagnold scaling for 
the shear stress S as 
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S = A(v)mo- 2 



- d i\ 



(1) 



where m is the grain mass, a is the grain diameter, and 
d is the spatial dimension, with A being a dimension- 
less coefficient that depends on v. Obviously this scaling 
should have broad range of validity for a simple shear 
flow of cohesionless hard grains; it should hold until ei- 
ther the system becomes so dense that the elasticity of 
particles comes into the problem or the shear banding 
destabilizes the uniform shear flow. In the case of gravi- 
tational flow down a slope, the shear is brought about by 
the gravity and the gravitational acceleration g brings 
another time scale into the problem, but the Bagnold 
scaling is expected to hold in the denser region where 
the effect of gravity on particle orbits between collisions 
is not significant. Actually the Bagnold scaling is ob- 
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FIG. 1: (color online) The y dependences of the packing frac- 
tion v (a), the granular temperature T (lines) and the rota- 
tional temperature T (marks) (b), and d y q/(S A /) (c) for var- 
ious inclination angles 8. The inset in (a) shows a schematic 
diagram of the system with a coordinate. For most of the 
data, the bottom boundary is BC1 (see text) and the total 
depth H is 50, but the data with BC2 and H — 100 are 
also given for 9 = 20°. In (b), T are divided by the factors 
0.456 (8 = 20°), 0.513 (8 = 21°), 0.571 (8 = 22°), and 0.696 
(8 — 23°) to show T oc T in the bulk region. The inset in (c) 
shows that T/7 2 is roughly constant in the bulk (see text). 



served in experiments 0, El and simulations in the 
slope flows quite well. 

Recently, Silbert et al. have performed large scale 
molecular dynamics simulations on dense slope flows 
and they found an interesting fact that the grain density 
is almost constant and independent of the depth except 
for the boundary layers near the bottom and the sur- 
face. This is quite intriguing because the density is not 
nearly the closed packed density; it depends upon the in- 
clination angle 9 but neither upon the total depth of the 
flow H 4] nor the roughness of the slope Somehow, 
the system adjusts its temperature to keep the density 
constant along the depth direction. 

This is, however, not difficult to understand if one ex- 
tends the Bagnold scaling to the pressure; under the same 
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condition, the pressure, or the normal stress N, should 
have the same form 

JV = B(iy)ma 2 - d -y 2 (2) 

as S with another dimensionless coefficient B, thus the 
ratio of S/N depends upon the packing fraction v, but 
not on the shear rate 7. In the gravitational slope flow, 
the force balance gives S/N — tan 6, thus we have 

A{v)/B{v) = tBXi6, (3) 

which shows the packing fraction is determined by the 
inclination 9 but does not depend upon the depth H. 
This argument suggests the existence of the bulk region 
with the density plateau in the gravitational flow is very 
general and independent of detailed properties of grains. 

In order to determine how the packing fraction v de- 
pends on 9, we need a theory that gives constitutive re- 
lations. This has been done by Louge |(| using a kinetic 
theory for inelastic hard spheres 0. It is disappointing, 
however, to find that the kinetic theory fails to give cor- 
rect density profiles; two branches of solution for v were 
found, but one gives too small v and the other gives oppo- 
site 9 dependence of v to the one observed in simulations, 
which implies the branch is a dynamically unstable one. 
This is a little puzzling situation because the kinetic the- 
ory has been shown to hold in the case of sheared flow in 
the similar density regime @ . 

In this paper, we present detailed analysis of our sim- 
ulations on the bulk region of two-dimensional gravita- 
tional flow, assuming the framework of kinetic theory. 
In contrast to previous works, where the overall profiles 
from hydrodynamic models were discussed Q, we exam- 
ine each constitutive relations separately using data in 
the bulk to avoid the uncertainty in a boundary condi- 
tion for hydrodynamic equations. 

First, we show how the bulk behavior is understood 
within the framework of kinetic theory. In the granu- 
lar kinetic theory, the kinetic temperature T = m < 
(c — v) 2 > jd is treated as a separate variable, which in- 
troduces an additional time scale. Here, c is the particle 
velocity, < • • • > represents average over the microscopic 
scales, and v =< c >. The shea r stre ss is given by the 
momentum flux; this is ml(y)jn\/T '/m in the elementary 
transport theory, where n is the number density and t{y ) 
represents the mean free path. More generally, 

S = / 2 (,)m 1 /V 1 -¥/ 2 i, (4) 

with a dimensionless function /2(f), which depends on 
v and other material parameters such as a restitution 
coefficient. Similarly we have for the normal stress N, 
the energy dissipation T, and the heat flux q 

N = h{v)o- d T, (5) 
T = h(v)m-^a- d -^l\ (6) 
q = -hMm-^a^T^dyT, (7) 



with d y = d/dy. The forms J3J to J7J) represent the 
quite general framework of kinetic theory [lj, although 
functional forms fi{v) vary depending upon level of ap- 
proximation. 

These expressions should be compatible with the Bag- 
nold scaling when the only relevant time scale is 7 . 
Equation Q of S indicates that T cx •j 2 in order that 
the Bagnold scaling should hold. In the kinetic the- 
ory, T is determined by the energy balance equation 

-d v q + 57 - r = (8) 

in the steady flow. When the divergence of the heat flux 
—d y q is zero as in the case of the uniform shear flow due 
to the symmetry in the y-direction, T is determined by 
the local balance between the viscous heating Sj and the 
energy loss T; then the time scale that determines T is 
the shear rate only. Sj — T with Eqs. 10} and © gives 

T=[f 2 (v)/f 3 (v)]ma 2 j 2 , (9) 

thus the Bagnold scaling holds. In the slope flow, the 
divergence of the heat flux is not necessarily zero, but 
it turns out to be small compared with the other terms. 
Therefore, from Eqs. QJ, JSJ, JSJ, and S/N = tan 6*, we 
have 

tan0= V/2M/3M//1M, (10) 

which gives v as a function of 9 if we know ji(y). 

In the following, the above analysis is examined in de- 
tail in comparison with the data of our two-dimensional 
simulations on the soft sphere model with the disk mass 
m, the diameter cr, and the moment of inertia I = 
mcr 2 /10 as in Ref. Q. The particles stiffness is taken 
to be in the region where the flow behavior is already 
in the hard sphere limit |10| . which allows us to em- 
ploy the constitutive relations for hard disks in the ki- 
netic theory analysis in the following. The linear spring- 
dashpot model and the Coulomb friction with the coef- 
ficient fi = 0.5 are employed, and the periodic boundary 
condition is imposed along the flow direction. The bot- 
tom boundary is made rough by attaching disks of di- 
ameter 2(T, which we refer to as BC1: See Ref.JH for de- 
tailed descriptions of the model (Our model corresponds 
to the model "L2"). We confirmed that our data agree 
with theirs in the bulk, although our slope length (20ct) 
is shorter than theirs (IOOct). We show only the data for 
9 > 20° with H — 50, which is well above the stopping 
angle 9 B t op (9 st0 p ~ 18° the boundary effects become 
significant for 9 closer to s t O p@]- The boundary effects 
are examined by simulating with a slope covered with 
disks of diameter a (BC2). 

To compare our data with the kinetic theory, we use 
the normal restitution coefficient e p = 0.92 and the tan- 
gential restitution coefficient j3 = 1, although the tan- 
gential restitution coefficient in the simulation is not con- 
stant because of sliding collisions with the Coulomb fric- 
tion The Coulomb friction is important in simula- 
tion, but no kinetic theories have been worked out yet 
with it in two-dimension 1121. 
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TABLE I: The constitutive relations from kinetic theory in Ref. |l3|. with parameters re = 4// (ma 2 ), a = re(l + /3)/[2(l + re)], 
r = (1 + e p )/2, and C = [-(1 - f/T)a 2 + (5 - 8r)a + 2(5 - 3r)| /2. g (u) is the radial distribution function. 
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Figure H shows the y-dependence of v (a), T (b), and 
d y q/(S A /) (c) for various inclination angles 8; most of the 
data are for the system with the depth H = 50 and BC1, 
but the data for BC2 and those for H = 100 are also 
shown for 8 = 20° for comparison. The data are given 
in the unit syst em where the length cr, the mass to, and 
the time \fa~j~g are unities. One can see that the packing 
fraction in the bulk does not depend on the depth, and 
the effects of the boundary condition are confined within 
the boundary layer and the bulk properties are indepen- 
dent. Figure ^c) shows that d y q is much smaller than 
Sj in the bulk, which is consistent with our argument to 
derive Eq. ©; the heat flux q is estimated by the consti- 
tutive relation in Ref.^3). The plots of T/7 2 in the inset 
shows that Eq. © holds approximately. 

We compare our data with the constitutive relations 
derived by Jenkins and Richman |l3j for two-dimensional 
inelastic hard disks. The functions /i(f), f?,{v), and 
/3(f) in the steady flow are given in Table D We adopt 
the radial distribution function go(v) from Ref. Q 



9o(v) = 9c(i 



9f{v)-9c{v) 
1 + exp(-(^ - v )/mo) ' 



(11) 



where g c (u) = (1 - 7i//16)(l - v)~ 2 and g f {u) = [(1 + 
e p )v(y/v c /v - I)]" 1 with v c = 0.82, u = 0.7006, and 
too = 0.0111. 

In /2(f) and fz{v), the rotational temperature T = 
I < (w — lo) 2 > appears as T/T, where w is the par- 
ticle angular velocity and u =< w >: In the kinetic 
theory, ui is simply assumed to be (V x v) z /2 13], which 
holds except for the region near the bottom boundary 
[l5|]. T/T becomes constant in the kinetic theory |l3j| : 
the value should be 1 for our parameters, but is not in the 
simulations. This should be mainly due to the Coulomb 
friction, which has strong effects on particle rotations. In 
Fig.QJb), T's are plotted by marks along with T's (lines): 
T's are divided by factors that give the best fits of T's 
with T's. We see T's fit with T's in the bulk region with 
the factor around 0.5, but the ratio depends on 8. In the 
following, we try both T/T — 1 and the values obtained 
from the simulations for T/T in f2{v) and fz(v). 

First, we examine Eqs. an d ©• Figure |21 shows 
the simulation data of N/T (a) and S/{jVT) (b) against 
v by marks. The data from the bottom layers (y < 10) 
are distinguished by filled marks because they follow a 
different trend. The data outside the bottom layers (y > 
10, open marks) for various 8 collapse onto a single line. 




FIG. 2: (color online) N/T (a) and S/{VTj) (b) vs. v 
for various 9. The open and filled marks represent the data 
outside and within the bottom boundary, respectively. The 
lines show fi(u) (a) and fiiy) (b). 
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FIG. 3: (color online) S^/T z/2 vs. v (marks) with f 3 (v) 
(lines) (a), and tan# vs. the bulk density (marks) with 
V /2M/3M//1M (lines) (b) for various 6. The lines rep- 
resent the plots with various values of T/T. 



This clearly shows that the expressions (0} and © in 
the kinetic theory are valid outside the bottom layers. 
The data from the bottom layers show some scatter and 
different tendency between BC1 and BC2. 

/i(z/) and f 2 [v) in Table[Qwith T/T = 1 are shown by 
the solid lines in Fig.^a) and (b), respectively, and they 
agree with the data. /2(f) depends on T/T only weakly 
and the difference turned out to be negligibly small in 
the range 0.5 < T/T < 1. 

Now, we examine fz(v) in Eq. ijfJJl. In Fig.|^a), we plot 
S7/T 3 / 2 against v from the data; this quantity should 
give fz(v) from Eqs. ijjjj and (0J. Only the data from the 
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bulk (15 < y < 35) are plotted because Eq. © is valid 
only in the bulk as we have already seen in Fig.^c). The 
lines show fs(y) from kinetic theory with various T/T. 
/3(f) depends on T/T, but the data agrees reasonably 
well with fa(v) when T/T ~ 0.5. Note, however, that 
the singularity at v = v c is weaker in the simulation data 
than in /3(f). 

This difference in /3(f) is significant when we see them 
in the bulk density. In Fig. EOT)), the bulk density v is 
plotted against the inclination angle 9 for the simula- 
tion data and for the kinetic theory; for the latter, we 
plot yj f2{v)}z{v) //1(f), which should give tan# from 
Eq. (|10|l . The bulk density decreases as 9 increases in the 
simulation, but yj (v)fs (y) / fi (y) shows opposite ten- 
dency; the density increases as 9 increases. This discrep- 
ancy comes mainly from the discrepancy in /3(f), more 
specifically from the fact that the data shows a weaker 
divergence in fi{v), while the kinetic theory assumes the 
same singularity in all of f\ (y) , $1 (y) , and fs (v) near the 
random closed packing. 

Some parts of the discrepancy might originate from the 
Coulomb friction, because it is included in the simulation 
and should have some effects on energy dissipation, but 
not taken into account in the existing two-dimensional 
theories. The existing three-dimensional theory [?J, how- 
ever, suggests that the way it changes fo{v) is just to 
modify the coefficient of v 2 gQ(y) as long as the level of 
approximation remains the same. Such a change is not 
enough to make the singularity in fo{v) weaker. 

It is a bit puzzling to find clear deviation from the ki- 
netic theory in the energy dissipation while the stresses 
agree quite well. A possible origin of this is the velocity 
correlation induced by the inelasticity, which could vio- 
late the molecular chaos assumption in the kinetic theory; 
The decrease of the relative velocity tends to reduce the 
energy loss per collision. This effect has been noticed 



in some granular gas simulations, where the energy loss 
rate is found to be more sensitive to the velocity corre- 
lation than stresses 16]. Careful analysis of the velocity 
correlation in dense flow is awaited. 

Before concluding, let us make some comments on the 
Pouliquen's flow rule Q: The flow velocity at the sur- 
face scales u{H)/ y/gH = bH / H stop (9) with H stop (9) be- 
ing the depth of the flow below which the flow stops for 
a given inclination angle 9, and a numerical constant b 
abound 0.136. Erta§ and Halsey ^} argued that the 
appearance of H stop (9) in the expression of flow velocity 
for the depth H, which can be much larger than H s t op (9) , 
implies that the rheology of the dense gravitational flow 
is not local, and have proposed the eddy mechanism. If 
the flow is controlled by a non-local mechanism, there 
is no way that the kinetic theory holds. The Pouliquen 
flow rule, however, does not necessarily means a non-local 
mechanism but it simply means the stopping depth is de- 
termined by some aspects of the flowing rheology. We do 
not know yet how it is determined, but the present re- 
sults suggest that the kinetic theory may well be a good 
starting point to describe the flow. 

In summary, by careful analysis of simulation data, 
we have demonstrated that the rheology of gravitational 
dense granular flow can be described within the frame- 
work of kinetic theory. Especially, the constitutive re- 
lations based on the kinetic theory have been shown to 
agree quantitatively with the simulations, but there is a 
slight discrepancy in the energy dissipation. Due to this 
discrepancy, the kinetic theory fails to give a correct de- 
scription of the density plateau of the granular flow down 
a slope. 
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